#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <opencv2/opencv.hpp>
//#include <Eigen/Eigen>
#include "fstream"
#include "dragon.h"
const float PI = 3.1415;
const int samplesize = 50000;
float getR(const float& r, const float& z_r, const float& A, const float& alpha2,const float& sigma_tr)
{
    float d_r = sqrt(r*r + z_r*z_r),z_v = z_r*(1+4*A/3), d_v = sqrt(r*r +z_v*z_v);
    float R;
    R = alpha2/(4*PI)*(z_r*(sigma_tr+1/d_r)*exp(-sigma_tr*d_r)/d_r/d_r +
                       z_v*(sigma_tr+1/d_v)*exp(-sigma_tr*d_v)/d_v/d_v);
    return R;
}
void getAllR(cv::Mat& R)
{
    float yita = 1.5;
    int samplesize = 50000;
    float sigma_s = 0.75, sigma_a = 0.02;
    float g=0.25;//散射角平均cos
    float sigma_s2 = (1-g)*sigma_s, sigma_t2 = sigma_a+sigma_s2,
            alpha2 = sigma_s2/sigma_t2, sigma_tr = sqrt(3*sigma_a*sigma_t2), z_r = 1/sigma_t2,
            F_dr = -1.440/yita*yita+0.710/yita+0.668+0.0636*yita, A = (1+F_dr)*(1-F_dr);
    for(int i=0;i<samplesize;i++){
        for(int j=0;j<samplesize;j++){
            float cur_r = sqrt(pow(vertices[i][0]-vertices[j][0],2)+
                    pow(vertices[i][1]-vertices[j][1],2)+
                    pow(vertices[i][2]-vertices[j][2],2));
            R.at<float>(i,j) = getR(cur_r,z_r,A,alpha2,sigma_tr);
        }
    }
}

int main()
{
    cv::Mat R(samplesize, samplesize, CV_32FC1);
    getAllR(R);
    std::cout<<R<<std::endl;
    std::ofstream out("R.txt",std::ios::binary);
    out.write((char*)R.data, R.rows * R.cols * sizeof(float));
    std::cout<<"output to file is finished"<<std::endl;
}
